%% FanChart
% This function generates a fan chart. A fan chart aims to visualize the
% uncertainty that surrounds a sequence of point fore{n}casts.
%
% Required input:
% -- a HxD matrix 'fore{n}': the density fore{n}casts D along the horizon 1:H
% -- a Tx1 vector 'hist{n}': a vector of hist{n}orical data
%
% Optional input:
% -- a string 'datefinal': referring to the end of the fore{n}cast horizon
% -- a scalar 'datestep': if e.g. set to 3, then labels quarterly 
% Varargin to use boxplots or not

function savename = FanChartmulti(foresets,legends,dimension,cfg,r,varargin)
global globals

if nargin == 5
    varargin{1} = true;
    varargin{2} = true;    
    varargin{3} = [0 1]/5;
elseif nargin == 6
    varargin{2} = true;
    varargin{3} = [0 1]/5;
elseif nargin == 7
    varargin{3} = [0 1]/5;
end

%% Figure
if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
globals.fid = figure(4);set(4,'Units','centimeter')
matrixsize = size(fore);

%% Boxplot

if varargin{1} && varargin{2}
    barhandle = subplot('position',[0.75 0.15 0.15 0.75]);
else
    barhandle = gca;
end
if varargin{2} 
configs.colors = hsl2rgb([((0:matrixsize(dimension)-1)*varargin{3}(2))+varargin{3}(1); ones([1,matrixsize(dimension)]); ones([1,matrixsize(dimension)])*0.6]');
aboxplot(permute(fore(:,1,end,:),[1 4 2 3]),'WidthS',1,'Labels',legends); m=0;
for n = [1:25]; child = get(gca,'Children'); child2=get(child(end),'Children'); try set(child2(n),'FaceColor',configs.colors(end-m,:)); m=m+1; catch; end; end
set(barhandle,'YColor',[1 1 1],'XColor',[1 1 1],'Ylim',[0 1].*get(barhandle,'Ylim'));
boxploth       = gca;
end

if varargin{1} && varargin{2}
fanh=subplot('position',[0.15 0.15 0.60 0.75]);
end

%% Magic
if varargin{2}
for n = matrixsize(dimension):-1:1
%% Set relevant quantiles (in percent)
% quantiles=[5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95];
quantiles  = linspace(25,75,19);
num_quant=size(quantiles,2);

mat_quant=NaN(matrixsize(3),num_quant);
for h=1:matrixsize(3) % loop over horizon
    for q=1:num_quant % loop over quantiles
        if dimension == 1
            mat_quant(h,q)=quantile(fore(n,:,h,:),quantiles(1,q)/100);
        elseif dimension == 4
            mat_quant(h,q)=quantile(fore(:,:,h,n),quantiles(1,q)/100);
        end
    end
end

colors  = hsl2rgb([ones([1,10])*((n-1)*varargin{3}(2)+varargin{3}(1)); ones([1,10]); linspace(0.6,0.35,9) 0.3]');
alpha   = .6;

hold on
for m = 1:9
    h = fill([1:matrixsize(3) flipdim(1:matrixsize(3),2)],...
        [mat_quant(:,m)' flipdim(mat_quant(:,m+1)',2) ],colors(m,:));
    set(h,'EdgeColor',colors(m,:),'EdgeAlpha',0,'FaceAlpha',alpha);
    
    h = fill([1:matrixsize(3) flipdim(1:matrixsize(3),2)],...
        [mat_quant(:,10*2-m)' flipdim(mat_quant(:,10*2-m-1)',2) ],colors(m,:));
    set(h,'EdgeColor',colors(m,:),'EdgeAlpha',0,'FaceAlpha',alpha);
end
mainh(n)=h;
m=10; h=fill([1:matrixsize(3) flipdim(1:matrixsize(3),2)],[mat_quant(:,m)' flipdim(mat_quant(:,m)',2)],colors(m,:));
set(h,'EdgeColor',colors(m,:),'FaceAlpha',0,'LineWidth',1.5);

end


%% Title+labels
title(['TravelDemand ' cfg.allregions{r} ', (n=' num2str(size(fore,1)) ')']);
if size(fore,3)==80
    PlotRange = 1:80;
else
    PlotRange = 35:130;
end
set(fanh,'YLim',get(barhandle,'YLim'),'XTick',cfg.Time(round(linspace(1,length(PlotRange),5))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(1,length(PlotRange),5))))));

ylabel('TravelDemand total, thousand pkm/yr');
xlabel('Year');
configs.ShowLegend=false;
if configs.ShowLegend
h=legend(mainh,legends,'Location','NorthWest'); 
end
pubgraph(globals.fid,10,2,'white');
end

if varargin{1}
set(boxploth,'YLim', get(gca,'YLim' ));
box(boxploth,'off');
end

savename = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}];

end

function rgb=hsl2rgb(hsl)

%Converts Hue-Saturation-Luminance Color value to Red-Green-Blue Color value
%
%Usage
%       RGB = hsl2rgb(HSL)
%
%   converts HSL, a M X 3 color matrix with values between 0 and 1
%   into RGB, a M X 3 color matrix with values between 0 and 1
%
%See also rgb2hsl, rgb2hsv, hsv2rgb

%Suresh E Joel, April 26,2003

if nargin<1,
    error('Too few arguements for hsl2rgb');
    return;
elseif nargin>1,
    error('Too many arguements for hsl2rgb');
    return;
end;

if max(max(hsl))>1 | min(min(hsl))<0,
    error('HSL values have to be between 0 and 1');
    return;
end;

for i=1:size(hsl,1),
    if hsl(i,2)==0,%when sat is 0
        rgb(i,1:3)=hsl(i,3);% all values are same as luminance
    end;
    if hsl(i,3)<0.5,
        temp2=hsl(i,3)*(1+hsl(i,2));
    else
        temp2=hsl(i,3)+hsl(i,2)-hsl(i,3)*hsl(i,2);
    end;
    temp1=2*hsl(i,3)-temp2;
    temp3(1)=hsl(i,1)+1/3;
    temp3(2)=hsl(i,1);
    temp3(3)=hsl(i,1)-1/3;
    for j=1:3,
        if temp3(j)>1, 
            temp3(j)=temp3(j)-1; 
        elseif temp3(j)<0, 
            temp3(j)=temp3(j)+1; 
        end;
        if 6*temp3(j)<1,
            rgb(i,j)=temp1+(temp2-temp1)*6*temp3(j);
        elseif 2*temp3(j)<1,
            rgb(i,j)=temp2;
        elseif 3*temp3(j)<2,
            rgb(i,j)=temp1+(temp2-temp1)*(2/3-temp3(j))*6;
        else
            rgb(i,j)=temp1;
        end;
    end;
end;

rgb=round(rgb.*100000)./100000; %Sometimes the result is 1+eps instead of 1 or 0-eps instead of 0 ... so to get rid of this I am rounding to 5 decimal places)
end

% 
% Copyright (C) 2011-2012 Alex Bikfalvi
%

function aboxplot(X,varargin)

% Parameters
widthl = 0.7;
widths = 0.8;
widthe = 0.4;
outMarker = '.';
outMarkerSize = 3;
outMarkerEdgeColor = [0.6 0.6 0.6];
outMarkerFaceColor = [0.6 0.6 0.6];
alpha = 0.05;
cmap = [];
colorrev = 0;
colorgrd = 'blue_down';

% Get the number or data matrices
if iscell(X)
    d = length(X);
else
    % If data is a matrix extend to a 3D array
    if 2 == ndims(X)
        X = reshape(X, [1,size(X)]);
    end
    d = size(X,1);
end;

% Get the data size
if iscell(X)
    n = size(X{1},2);
else
    n = size(X,3);
end

% Set the labels
labels = cell(n,1);
for i=1:n
    labels{i} = num2str(i);
end

% Optional arguments
optargin = size(varargin,2);

i = 1;
while i <= optargin
    switch lower(varargin{i})
        case 'labels'
            labels = varargin{i+1};
        case 'colormap'
            cmap = varargin{i+1};
        case 'colorgrad'
            colorgrd = varargin{i+1};
        case 'colorrev'
            colorrev = varargin{i+1};
        case 'outliermarker'
            outMarker = varargin{i+1};
        case 'outliermarkersize'
            outMarkerSize = varargin{i+1};
        case 'outliermarkeredgecolor'
            outMarkerEdgeColor = varargin{i+1};
        case 'outliermarkerfacecolor'
            outMarkerFaceColor = varargin{i+1};
        case 'widthl'
            widthl = varargin{i+1};
        case 'widths'
            widths = varargin{i+1};
        case 'widthe'
            widthe = varargin{i+1};
    end
    i = i + 2;
end

% Colors
colors = cell(d,n);

if colorrev
    %  Set colormap
    if isempty(cmap)
        cmap = colorgrad(n,colorgrd);
    end
    if size(cmap,1) ~= n
        error('The number of colors in the colormap must equal n.');
    end
    for j=1:d
        for i=1:n
            colors{j,i} = cmap(i,:);
        end
    end
else
    %  Set colormap
    if isempty(cmap)
        cmap = colorgrad(d,colorgrd);
    end
    if size(cmap,1) ~= d
        error('The number of colors in the colormap must equal n.');
    end
    for j=1:d
        for i=1:n
            colors{j,i} = cmap(j,:);
        end
    end
end

xlim([0.5 n+0.5]);

hgg = zeros(d,1);

for j=1:d
    % Get the j matrix
    if iscell(X)
        Y = X{j};
    else
        Y = squeeze(X(j,:,:));
    end
    
    % Create a hggroup for each data set
    hgg(j) = hggroup();
    set(get(get(hgg(j),'Annotation'),'LegendInformation'),'IconDisplayStyle','on');
    legendinfo(hgg(j),'patch',...
        'LineWidth',0.5,...
        'EdgeColor','k',...
        'FaceColor',colors{j,1},...
        'LineStyle','-',...
        'XData',[0 0 1 1 0],...
        'YData',[0 1 1 0 0]);
    
    for i=1:n
    
        % Calculate the mean and confidence intervals
        [q1 q2 q3 fu fl ou ol] = quartile(Y(:,i));
        u = nanmean(Y(:,i));

        % large interval  [i - widthl/2 i + widthl/2] delta = widthl
        % medium interval start: i - widthl/2 + (j-1) * widthl / d
        % medium interval end: i - widthl/2 + j * widthl / d
        % medium interval width: widthl / d
        % medium interval middle: i-widthl/2+(2*j-1)*widthl/(2*d)
        % small interval width: widths*widthl/d
        % small interval start: i-widthl/2+(2*j-1-widths)*widthl/(2*d)
  
        % Plot outliers
        hold on;
        plot((i-widthl/2+(2*j-1)*widthl/(2*d)).*ones(size(ou)),ou,...
            'LineStyle','none',...
            'Marker',outMarker,...
            'MarkerSize',outMarkerSize,...
            'MarkerEdgeColor',outMarkerEdgeColor,...
            'MarkerFaceColor',outMarkerFaceColor,...
            'HitTest','off',...
            'Parent',hgg(j));
        plot((i-widthl/2+(2*j-1)*widthl/(2*d)).*ones(size(ol)),ol,...
            'LineStyle','none',...
            'Marker',outMarker,...
            'MarkerSize',outMarkerSize,...
            'MarkerEdgeColor',outMarkerEdgeColor,...
            'MarkerFaceColor',outMarkerFaceColor,...
            'HitTest','off',...
            'Parent',hgg(j));
        hold off;
        
        % Plot fence
        line([i-widthl/2+(2*j-1)*widthl/(2*d) i-widthl/2+(2*j-1)*widthl/(2*d)],[fu fl],...
            'Color','k','LineStyle',':','HitTest','off','Parent',hgg(j));
        line([i-widthl/2+(2*j-1-widthe)*widthl/(2*d) i-widthl/2+(2*j-1+widthe)*widthl/(2*d)],[fu fu],...
            'Color','k','HitTest','off','Parent',hgg(j));
        line([i-widthl/2+(2*j-1-widthe)*widthl/(2*d) i-widthl/2+(2*j-1+widthe)*widthl/(2*d)],[fl fl],...
            'Color','k','HitTest','off','Parent',hgg(j));
        
        % Plot quantile
        if q3 > q1
            rectangle('Position',[i-widthl/2+(2*j-1-widths)*widthl/(2*d) q1 widths*widthl/d q3-q1],...
                'EdgeColor','k','FaceColor',colors{j,i},'HitTest','off','Parent',hgg(j));
        end
        
        % Plot median
        line([i-widthl/2+(2*j-1-widths)*widthl/(2*d) i-widthl/2+(2*j-1+widths)*widthl/(2*d)],[q2 q2],...
            'Color','k','LineWidth',1,'HitTest','off','Parent',hgg(j));
        
        % Plot mean
        hold on;
        plot(i-widthl/2+(2*j-1)*widthl/(2*d), u,...
            'LineStyle','none',...
            'Marker','o',...
            'MarkerEdgeColor','k',...
            'MarkerFaceColor',colors{j,i},...
            'HitTest','off','Parent',hgg(j));
        hold off;
    end
end

box on;

set(gca,'XTick',1:n);
set(gca,'XTickLabel',labels);

end

function c = colorgrad(varargin)

n = 16;
t = 'blue_down';

switch length(varargin)
    case 1
        n = varargin{1};
    case 2
        n = varargin{1};
        t = varargin{2};
end

switch lower(t)
    case 'blue_up'
        c = cat(2,linspace(0,0.6,n)',linspace(0.2,0.8,n)',linspace(0.6,1,n)');
    case 'blue_down'
        c = cat(2,linspace(0.6,0,n)',linspace(0.8,0.2,n)',linspace(1,0.6,n)');
    case 'orange_up'
        c = cat(2,linspace(1,248/255,n)',linspace(0.6,224/255,n)',linspace(0,124/255,n)');
    case 'orange_down'
        c = cat(2,linspace(248/255,1,n)',linspace(224/255,0.6,n)',linspace(124/255,0,n)');
    case 'green_up'
        c = cat(2,linspace(0.2,0.6,n)',linspace(0.6,1,n)',linspace(0.2,0.6,n)');
    case 'green_down'
        c = cat(2,linspace(0.6,0.2,n)',linspace(1,0.6,n)',linspace(0.6,0.2,n)');
    case 'red_up'
        c = cat(2,linspace(.8,1,n)',linspace(.2,.6,n)',linspace(.2,.6,n)');
    case 'red_down'
        c = cat(2,linspace(1,.8,n)',linspace(.6,.2,n)',linspace(.6,.2,n)');
    otherwise
        error('No such color gradient.');
end

end


function [q1 q2 q3 fu fl ou ol] = quartile(x)

% rank the data
y = sort(x);

% compute 50th percentile (second quartile)
q2 = quantile(y,0.5);

% compute 25th percentile (first quartile)
q1 = quantile(y,0.25);

% compute 75th percentile (third quartile)
q3 = quantile(y,0.75);

% rank the data
y = sort(x);

% compute 50th percentile (second quartile)
q2 = nanmedian(y);

% compute 25th percentile (first quartile)
q1 = nanmedian(y(y<=q2));

% compute 75th percentile (third quartile)
q3 = nanmedian(y(y>=q2));

% compute Interquartile Range (IQR)
IQR = q3-q1;

fl = min(y(y>=q1-1.5*IQR));
fu = max(y(y<=q3+1.5*IQR));

ol = y(y<q1-1.5*IQR);
ou = y(y>q3+1.5*IQR);

end